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Abstract. We perform non-equilibrium simulations to study heat conduction in two-dimensional 
strongly coupled dusty plasmas. Temperature gradients are established by heating one part of the 
otherwise equilibrium system to a higher temperature. Heat conductivity is measured directly from the 
stationary temperature profile and heat flux. Particular attention is paid to the influence of damping 
effect on the heat conduction. It is found that the heat conductivity increases with the decrease of the 
damping rate, while its magnitude agrees with previous experimental measurement. 

PACS numbers: 52.25.Fi, 52.27.Gr, 52.27.Lw 

Introduction: Recently, experiments had been carried out to study thermal conduction 
in two-dimensional (2D) strongly coupled dusty plasmas (SCDPs) [H [2] in both 
crystalline and solid/liquid mixture states, and a thermal conductivity, which is 
independent of temperature, was found. Although both these experiments were aimed 
at studying the heat conduction at an atomic (molecular) level, neither of them showed 
many details of microscopic processes during the heat transfer. Therefore, we conduct 
here non-equilibrium simulation by using Brownian Dynamics method [3] to study heat 
transfer in 2D SCDPs in more details, serving as a supplement to real experiment. 

Numerical simulation: iV = 10000 particles are simulated in a rectangular area with 
periodical boundary condition in y direction and confining boundary condition in x 
direction. (More details of simulation and algorithm may be found in Ref. [3].) Particles 
interact with each other via pairwise Yukawa potential: <f>(r) = {Q 2 /r) exp (— r/Xp), 
with Q, r and being the particle charge, interparticle-distance and screening length, 
respectively. The strong-coupling strength is given by T = Q 2 /(ak B T), and the screening 
parameter by k = a/Xo, where a = (irn)^ 1 ^ 2 is the 2D Wigner-Seitz radius with n being 
the areal number density and k B T the system temperature. In addition, the damping 
coefficient 7 is needed to fully characterize the dynamics of the system. To simplify later 
discussion, we also introduce here the nominal plasma frequency co = [2<5 2 /(ma 3 )] 1 / 2 , 
where m is the mass of a particle. In the simulation, the screening parameter is kept 
constant at k = 1, as it is the most typical value found in experiment, while T and 7 
are varied to realize different equilibrium states and different damping rates. 

Our simulation is directly mimicking recent experiments [HE], and is different from 
the usual method of non-equilibrium simulation for heat conduction [I]. The system 
is firstly brought to an equilibrium with desired temperature (To) in either liquid or 
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Figure 1. Stationary temperature profiles for: (a) 7 = 0.05wo but different system states and 
temperature gradients and; (b) To = 0.36T* and T\ = 3.6T* but different damping rate. In both figures, 
symbols are simulation results while solid lines are fits according to Eq. ([1]). For (a), A = 0.35a;ofc_B, 
and for (b) different A values (given in the inserted plot) are obtained. 

solid state. The melting point for k = 1 is at T* fa 180 [5], and we'll denote the 
corresponding temperature as T*. Then the right half of the system (x > 0), is heated 
to a higher temperature (Ti) by applying a Gaussian white noise with desired strength. 
The evolution of the temperature profile and also the heat flux are recorded. A steady 
state is approached after a substantially long period. 

The microscopic heat flux for ith particle is defined as: J«(t) = Vj-Ej + 
\ Ylj=i,j^i r ij(^ij ■ v ij) - r i( F ext ■ Vj) where E { = (l/2)(mvf + J2j=i,j^i<f>ij) + <f>ext is 
the particle energy. The total flux in a region is then a summation of the microscopic 
flux of all particles therein divided by its area A, i. e., J(t) = (1/A) J2ieA We 
are mainly interested in the x-component of the heat flux J x . The three terms on the 
right-hand-side of above equation correspond to respectively contributions from: (1) 
the particle migration, which is believed to be the main mechanism of heat transport in 
gas and is denoted as Jk x hereafter, (2) particle interactions, i. e., phonon scattering, 
which is dominant in solids and is denoted as J px hereafter, and (3) the external force. 
Since the external force acts on only a few rows of particles around the two confining 
boundaries, its direct influence on the heat flux is localized. Neglecting the external 
contribution will bring it back to the standard one (6j [7] . 

Analytical model: The heat transfer in our specific case can be described by Fourier's 
law: J = — AVT together with energy balance between heat conduction and energy 
dissipation due to damping: V(AVT) = 2 r yn(T — T )kB [2], where A is the heat 
conductivity. One has, 

T(x)-T = ^ T ^eV^ Ji -,(a;<0); T(x) - T x = ^-^TV^ 1 *, ( x > 0). (1) 

Expressions for heat flux may be obtained in a straightforward way, and we omit the 
results here. Since A is the only unknown parameter in Eq. ([T]), it may be measured 
by fitting the stationary temperature profile (STP) to Eq. ([T|). It should also be 
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mentioned that the Fourier's law could break down for low dimensional crystalline 
systems [TJ, largely due to a slow decay of equilibrium correlations of the heat current 
and a divergence of the finite-size conductivity. However, both of them may be avoided 
in dusty plasmas because of the finite damping effect. Therefore we skip this question 
at this moment, while interested readers may find more discussions in [7J. 

Results and discussions: Figure Q] (a) shows examples of STPs for different system 
states and temperature gradients with 7 = 0.05u;o, which is close to the experimental 
condition of [2J. Symbols are measurements from simulation, while solid lines are 
analytical fits according to Eq. ([TJ). These fits give a constant heat conductivity, 
A = (0.35 ± 0.05)u;o^B) or i n terms of thermal diffusivity Dj<: ~ 22mm 2 /s using 
parameters from Ref. [2j. This value is between the experimental measurement for 
crystalline state (30mm 2 /s) [1] and that for solid/liquid mixture phase (9mm 2 /s) [2]. 
Fits for high temperature (e. g. the two upper-most curves) suggest a slightly smaller 
A. Nevertheless, the value is in the range of the error bar for the present measurement. 

Figure [1] (b) shows STPs for different damping rate with other parameters fixed. 
Fits with Eq. (JTJ give a damping- dependent heat conductivity, as is shown in the 
inserted plot that A rises slightly with the decrease of 7. Note that this tendency 
is contradictory with that given by the analytical model in [8], which predicts an 
increase of A with increase of damping rate and was confirmed by their experiment 
[8]. However, this model is based on an empirical relation between diffusion and heat 
conduction coefficients obtained by fitting simulation results for three-dimensional (3D) 
simple liquid without damping, and their experiment was also performed in a 3D dusty 
plasma liquid [Bj. Whereas in our simulation we study 2D systems covering both liquid 
and solid states and/or with a liquid-solid mixture phase. So the discrepancy could 
have been caused by the different dimensionality and system states, as is known that 
transport processes depends much on these two factors. In our simulation, the damping 
effect is taken into account self-consistently and increase of A with decrease of 7 may be 
intuitively understood as follows. It has two effects on heat conduction: direct energy 
dissipation and indirect suppression of phonon propagation. The first one is only related 
to kinetic energy of the system, and had been explicitly taken into account in Eq. (CQ), 
whereas the second one affects the collective modes and is not included in Eq. (JTJ). 
Therefore decrease of 7 means less damping of phonon propagation, more efficient heat 
transfer through phonon scattering and consequently a higher heat conductivity. 

Figure [2] shows distributions of J x , J px and J^ x for different system states and 
temperature gradients, together with the analytical result derived from Eq. (CQ). Firstly, 
it may be seen that the kinetic part Jk x and the phonon part J px have different weights 
in different system states. As expected, J px is clearly dominant for solid state and low 
temperature liquid state, while Jk x dominates for high temperature liquid. The critical 
temperature where the two parts become equal is about 6T*. Secondly, the heat flux 
is not symmetrical about the heating interface. The decay of heat flux on the low- 
temperature side is slower and one needs a smaller slope, consequently a larger A to fit 
J x on this side, indicating a higher heat conductivity for lower temperature. Thirdly, the 
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Figure 2. Heat flux for different system states and different temperature gradients with 7 = 0.05loq. 
Symbols are simulation results, solid lines are fits according to Eq. {T]), using A = 0.35wofcs, an d dash 
lines in the two panels on the left are direct linear fits of the heat flux with A = 1.2o?oA>£. 



agreement between analytical results and simulation depends closely on system states. 
One generally observes a better agreement on the high temperature side and for higher 
temperature. These features suggest that A becomes temperature-dependent. 

Thus, we have measured the heat conductivity A of 2D SCDPs by analyzing both 
stationary temperature profile and heat flux in non-equilibrium simulations. It is found 
that A increases with the decrease of the damping rate. In addition, our results also 
suggest that A should be temperature-dependent. 
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